!------------------------------------------------------------------------------
! TJU/Department of Mechanics, Fluid Mechanics, Code START
!------------------------------------------------------------------------------
!> @Author: Jianxin Liu shookware@tju.edu.cn
!> @Date: 2024-08-15 22:58:10
!> @LastEditTime: 2024-08-15 22:58:11
!> @LastEditors: Jianxin Liu shookware@tju.edu.cn
!> @Description:
!> @FilePath: /START_NS_3D/src/reconstructor/module_reconstructor_base.f90
!> @Copyright (c) 2024 by Jianxin Liu email: shookware@tju.edu.cn, All Rights Reserved.
!------------------------------------------------------------------------------
module mod_reconstructor

   use petsc_wrapper
   use stdlib_optval
   use iso_fortran_env, only: dp => real64

   private
   public :: type_reconstructor_base
   public :: type_reconstructor_upwind_7order
   public :: type_reconstructor_linear_tvd
#include <start.h90>
   type, abstract :: type_reconstructor_base
      integer :: length
      integer :: stencil_width
      logical :: flg_bounds(2)
      integer :: ista, iend

   contains
      procedure(create_interface), deferred :: create
      procedure(reconstruct_interface), deferred :: reconstruct
      procedure(destroy_interface), deferred :: destroy

   end type type_reconstructor_base

   abstract interface

      subroutine create_interface(this, length, stencil_width, flg_bounds, limiter)
         import
         class(type_reconstructor_base), intent(inout) :: this
         integer, intent(in) :: length
         integer, intent(in) :: stencil_width
         logical, intent(in) :: flg_bounds(2)
         integer, intent(in), optional :: limiter
      end subroutine create_interface

      subroutine reconstruct_interface(this, f, h, s)
         import
         class(type_reconstructor_base), intent(in) :: this
         class(type_vector_base), intent(in) :: f
         class(type_vector_base), intent(inout) :: h
         character, intent(in), optional :: s
      end subroutine reconstruct_interface

      subroutine destroy_interface(this)
         import
         class(type_reconstructor_base), intent(inout) :: this
      end subroutine destroy_interface

   end interface

   type, extends(type_reconstructor_base), abstract :: type_reconstructor_upwind

   end type type_reconstructor_upwind

   type, extends(type_reconstructor_upwind) :: type_reconstructor_upwind_7order

   contains
      procedure :: create => CreateUpwind7Order
      procedure :: reconstruct => ReconstructUpwind7Order
      procedure :: destroy => DestroyUpwind7Order
   end type type_reconstructor_upwind_7order

   interface
      module subroutine CreateUpwind7Order(this, length, stencil_width, flg_bounds, limiter)
         class(type_reconstructor_upwind_7order), intent(inout) :: this
         integer, intent(in) :: length
         integer, intent(in) :: stencil_width
         logical, intent(in) :: flg_bounds(2)
         integer, intent(in), optional :: limiter
      end subroutine CreateUpwind7Order

      module subroutine ReconstructUpwind7Order(this, f, h, s)
         class(type_reconstructor_upwind_7order), intent(in) :: this
         class(type_vector_base), intent(in) :: f
         class(type_vector_base), intent(inout) :: h
         character, intent(in), optional :: s
      end subroutine ReconstructUpwind7Order

      module subroutine DestroyUpwind7Order(this)
         class(type_reconstructor_upwind_7order), intent(inout) :: this
      end subroutine DestroyUpwind7Order
   end interface

   type, extends(type_reconstructor_base), abstract :: type_reconstructor_tvd

      integer :: scheme_limiter
      type(type_idxset_general) :: is_p
      type(type_idxset_general) :: is_m
      procedure(limiter_interface), pointer :: limiter => null()

   contains
      procedure :: select_limiter => SelectLimiterTVD

   end type type_reconstructor_tvd

   abstract interface
      pure function limiter_interface(this, r) result(phi)
         import
         class(type_reconstructor_tvd), intent(in) :: this
         real(dp), intent(in) :: r
         real(dp) :: phi
      end function limiter_interface
   end interface

   type, extends(type_reconstructor_tvd) :: type_reconstructor_linear_tvd

   contains
      procedure :: create => CreateLinearTVD
      procedure :: reconstruct => ReconstructLinearTVD
      procedure :: reconstruct2 => ReconstructLinearTVD2
      procedure :: destroy => DestroyUpLinearTVD
   end type type_reconstructor_linear_tvd

   type, extends(type_reconstructor_tvd) :: type_reconstructor_muscl
      real(dp) :: kappa = 1.0d0/3.0d0
   contains
      procedure :: create => CreateMUSCL
      procedure :: reconstruct => ReconstructMUSCL
      procedure :: reconstruct2 => ReconstructMUSCL2
      procedure :: destroy => DestroyUpMUSCL

   end type type_reconstructor_muscl

   type, abstract, extends(type_reconstructor_base) :: type_reconstructor_weno

   contains

   end type type_reconstructor_weno

   type, extends(type_reconstructor_weno) :: type_reconstructor_weno_5order

   contains
      procedure :: create => CreateWENO5Order
      procedure :: reconstruct => ReconstructWENO5Order
      procedure :: reconstruct2 => ReconstructWENO5Order2
      procedure :: destroy => DestroyUpWENO5Order
   end type type_reconstructor_weno_5order

   #: set LIMITER_NAMES=["LinearTVD", "MUSCL", "WENO5Order"]
   #: set LIMITER_TYPES=["linear_tvd", "muscl", "weno_5order"]
   #: for limiter_name, limiter_type in list(zip(LIMITER_NAMES, LIMITER_TYPES))
   interface
      module subroutine Create${limiter_name}$ (this, length, stencil_width, flg_bounds, limiter)
         class(type_reconstructor_${limiter_type}$), intent(inout) :: this
         integer, intent(in) :: length
         integer, intent(in) :: stencil_width
         logical, intent(in) :: flg_bounds(2)
         integer, intent(in), optional :: limiter
      end subroutine Create${limiter_name}$

      module subroutine Reconstruct${limiter_name}$ (this, f, h, s)
         class(type_reconstructor_${limiter_type}$), intent(in) :: this
         class(type_vector_base), intent(in) :: f
         class(type_vector_base), intent(inout) :: h
         character, intent(in), optional :: s
      end subroutine Reconstruct${limiter_name}$

      module subroutine DestroyUp${limiter_name}$ (this)
         class(type_reconstructor_${limiter_type}$), intent(inout) :: this
      end subroutine DestroyUp${limiter_name}$

      module subroutine Reconstruct${limiter_name}$2(this, f, h, s)
         class(type_reconstructor_${limiter_type}$), intent(in) :: this
         class(type_self_vector), intent(in) :: f
         class(type_self_vector), intent(inout) :: h
         character, intent(in), optional :: s
      end subroutine Reconstruct${limiter_name}$2
   end interface
   #: endfor

contains

   subroutine SelectLimiterTVD(this, limiter)
      implicit none
      class(type_reconstructor_tvd), intent(inout) :: this
      integer, intent(in) :: limiter
      select case (limiter)
      case (LIMITER_MINMOD)
         this%limiter => minmod
      case (LIMITER_VAN_LEER)
         this%limiter => van_leer
      case (LIMITER_VAN_ALBADA)
         this%limiter => van_albada
      case (LIMITER_NONE)
         this%limiter => none_limiter
      case (LIMITER_SUPERBEE)
         this%limiter => superbee
      case (LIMITER_DUAL_MINMOD)
         this%limiter => dual_minmod
      case default
         call petsc_error_print("unknown limiter.")
      end select

   end subroutine SelectLimiterTVD

   pure function minmod(this, r) result(phi)
      implicit none
      class(type_reconstructor_tvd), intent(in) :: this
      real(dp), intent(in) :: r
      real(dp) :: phi

      phi = min(1.0_dp, max(r, 0.0_dp))

   end function minmod

   pure function van_leer(this, r) result(phi)
      implicit none
      class(type_reconstructor_tvd), intent(in) :: this
      real(dp), intent(in) :: r
      real(dp) :: phi

      phi = (r + abs(r))/(1.0_dp + r)

   end function van_leer

   pure function van_albada(this, r) result(phi)
      implicit none
      class(type_reconstructor_tvd), intent(in) :: this
      real(dp), intent(in) :: r
      real(dp) :: phi

      phi = (r + 1.0_dp)/(1.0d0/r + r)

   end function van_albada

   pure function none_limiter(this, r) result(phi)
      implicit none
      class(type_reconstructor_tvd), intent(in) :: this
      real(dp), intent(in) :: r
      real(dp) :: phi
      phi = 0.0_dp
   end function none_limiter

   pure function superbee(this, r) result(phi)
      implicit none
      class(type_reconstructor_tvd), intent(in) :: this
      real(dp), intent(in) :: r
      real(dp) :: phi

      phi = max(min(2.0_dp*r, 1.0_dp), min(r, 2.0_dp))

   end function superbee

   pure function dual_minmod(this, r) result(phi)
      implicit none
      class(type_reconstructor_tvd), intent(in) :: this
      real(dp), intent(in) :: r
      real(dp) :: phi

      if (r > 0.0_dp) then
         phi = min(2.0_dp*r, 2.0_dp, (1.0_dp + r)/2.0_dp)
      else
         phi = 0.0_dp
      end if

   end function dual_minmod

end module mod_reconstructor
